In [2]:
## Import Modules
from sklearn.feature_selection import RFECV
from sklearn.feature_selection import RFE
import statsmodels.api as sm
from sklearn.metrics import completeness_score, homogeneity_score
from sklearn.cluster import KMeans
from sklearn.model_selection import GridSearchCV
import math
from sklearn import feature_selection
import graphviz
from sklearn import metrics
from sklearn.tree import export_graphviz
from sklearn import cross_validation
from sklearn.metrics import classification_report
from sklearn.metrics import confusion_matrix
from sklearn import tree
import seaborn as sns
from sklearn import preprocessing
import numpy as np
import pandas as pd
import os
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
3#from sklearn import tree, naive_bayes
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
In [3]:
##get working directory
os.getcwd()
Out[3]:
In [4]:
os.chdir('C:\\Depaul_Win7\\DSC 478 Machine Learning\\Final Project\\health_adult')
In [5]:
##load data via read_csv
health_df = pd.read_csv('samadult.csv', sep=",")
health_df.head()
Out[5]:
In [6]:
#check types
health_df.dtypes
Out[6]:
In [7]:
columnNames = list(health_df.head(0))
print(columnNames)
In [8]:
#create reduced dataframe
In [9]:
#create reduced dataframe
health_red = health_df[['SEX', 'R_MARITL', 'MRACRPI2', 'REGION', 'PAR_STAT', 'DOINGLWA', 'SUPERVIS', 'WRKCATA', 'HOURPDA', 'PDSICKA', 'WRKLYR4', 'HYPEV', 'HYBPLEV', 'CHLEV', 'CHDEV', 'MIEV', 'STREV', 'COPDEV', 'AASMEV', 'ULCEV', 'CANEV', 'DBHVPAY', 'DBHVCLY', 'DBHVWLY', 'DBHVPAN', 'DBHVCLN', 'DBHVWLN', 'DIBREL', 'DIBEV1', 'DIBPRE2', 'EPILEP1', 'AHAYFYR', 'SINYR', 'CBRCHYR', 'KIDWKYR', 'LIVYR', 'JNTSYMP', 'ARTH1', 'PAINECK', 'PAINLB', 'PAINFACE', 'AMIGR', 'ACOLD2W', 'AINTIL2W', 'AHEARST1', 'AVISION', 'VIM_GLEV', 'VIM_MDEV', 'VIMGLASS', 'AVISACT', 'CHPAIN6M', 'AHSTATYR', 'FLA1AR', 'SPECEQ', 'ALC1YR', 'CIGAREV2', 'ECIGEV2', 'SMKSTAT2', 'APLKIND', 'AWORPAY', 'ADNLONG2', 'ASRGYR', 'AMDLONGR', 'HIT1A', 'HIT2A', 'HIT3A', 'HIT4A', 'FLUVACYR', 'LIVEV', 'ASICPUSE', 'ASIRETR', 'ASIMEDC', 'ASISTLV', 'ASICNHC', 'AWEBUSE', 'YTQU_YG1','ALCSTAT','YRSWRKPA', 'ASISLEEP', 'AHEIGHT', 'AWEIGHTP', 'BMI', 'BEDDAYR', 'CLCKTP', 'AGE_P', 'AHCNOYR2', 'LOCALL1B']]
In [10]:
#recode race
health_red['MRACRPI2'] = health_red['MRACRPI2'].replace([9,10,11,15,16,17], 4)
In [11]:
health_red['MRACRPI2'].value_counts().plot(kind='bar')
Out[11]:
In [12]:
health_red.columns.get_loc("R_MARITL")
Out[12]:
In [13]:
#recode marital status
health_red['R_MARITL'] = health_red['R_MARITL'].replace([2,3], 1) #married
In [14]:
health_red['R_MARITL'] = health_red['R_MARITL'].replace([4], 2) #widowed
In [15]:
health_red['R_MARITL'] = health_red['R_MARITL'].replace([5,6], 3) #divorced
In [16]:
health_red['R_MARITL'] = health_red['R_MARITL'].replace([7,8,9], 4) #never married
In [17]:
health_red['R_MARITL'].value_counts().plot(kind='bar')
Out[17]:
In [18]:
#recode missing values in part of dataframe
health_nan = health_red.iloc[:,0:76].replace([7,8,9], np.nan)
health_nan2 = health_red.iloc[:,76]
health_nan3 = pd.concat([health_nan,health_nan2],axis=1,join_axes=[health_nan.index])
health_nan3.head()
Out[18]:
In [19]:
health_nan3.shape
Out[19]:
In [20]:
#recode values in rest of dataframe to merge
health_rest = health_red.iloc[:,77:87].replace([96,97,98,99,996,997,998,999,9999], np.nan)
health_rest.head()
Out[20]:
In [21]:
health_rest.shape
Out[21]:
In [22]:
#merge dataframes with values 7,8,9,96,97,98,99,996,997,998,999,9999 recoded to NaN
mergedf=pd.concat([health_nan3,health_rest],axis=1,join_axes=[health_nan3.index])
mergedf.head()
Out[22]:
In [23]:
mergedf['LOCALL1B'].shape
Out[23]:
In [24]:
#test.isna().sum()
print(mergedf.isna().sum().to_string())
In [25]:
#remove NAs
health_na = mergedf.dropna()
In [26]:
#check number of NAs
health_na.isnull().sum().sum()
Out[26]:
In [27]:
health_na.shape
Out[27]:
In [28]:
health_na['LOCALL1B'].shape
Out[28]:
In [29]:
health_na['ASISLEEP'].value_counts().plot(kind='bar')
Out[29]:
In [30]:
health_na['ASISLEEP'].hist()
Out[30]:
In [31]:
health_na.describe()
Out[31]:
In [32]:
health_na[health_na.columns[0]].value_counts().plot(kind='bar')
Out[32]:
In [33]:
health_na[health_na.columns[1]].value_counts().plot(kind='bar')
Out[33]:
In [34]:
health_na[health_na.columns[2]].value_counts().plot(kind='bar')
Out[34]:
In [35]:
health_na[health_na.columns[3]].value_counts().plot(kind='bar')
Out[35]:
In [36]:
health_na[health_na.columns[4]].value_counts().plot(kind='bar')
Out[36]:
In [37]:
health_na[health_na.columns[5]].value_counts().plot(kind='bar')
Out[37]:
In [38]:
#plot all variables
# matplotlib.style.use('ggplot')
# health_na.plot.bar(subplots=True)
In [39]:
#save data frame without NAs to csv file
health_na.to_csv('health_na.csv', sep=",")
In [40]:
#create reduced categorical dataframe
health_cat = health_na[['SEX', 'R_MARITL', 'MRACRPI2', 'REGION', 'PAR_STAT', 'DOINGLWA', 'SUPERVIS', 'WRKCATA', 'HOURPDA', 'PDSICKA', 'WRKLYR4', 'HYPEV', 'HYBPLEV', 'CHLEV', 'CHDEV', 'MIEV', 'STREV', 'COPDEV', 'AASMEV', 'ULCEV', 'CANEV', 'DBHVPAY', 'DBHVCLY', 'DBHVWLY', 'DBHVPAN', 'DBHVCLN', 'DBHVWLN', 'DIBREL', 'DIBEV1', 'DIBPRE2', 'EPILEP1', 'AHAYFYR', 'SINYR', 'CBRCHYR', 'KIDWKYR', 'LIVYR', 'JNTSYMP', 'ARTH1', 'PAINECK', 'PAINLB', 'PAINFACE', 'AMIGR', 'ACOLD2W', 'AINTIL2W', 'AHEARST1', 'AVISION', 'VIM_GLEV', 'VIM_MDEV', 'VIMGLASS', 'AVISACT', 'CHPAIN6M', 'AHSTATYR', 'FLA1AR', 'SPECEQ', 'ALC1YR', 'CIGAREV2', 'ECIGEV2', 'SMKSTAT2', 'APLKIND', 'AWORPAY', 'ADNLONG2', 'ASRGYR', 'AMDLONGR', 'HIT1A', 'HIT2A', 'HIT3A', 'HIT4A', 'FLUVACYR', 'LIVEV', 'ASICPUSE', 'ASIRETR', 'ASIMEDC', 'ASISTLV', 'ASICNHC', 'AWEBUSE', 'YTQU_YG1', 'ALCSTAT']]
In [41]:
health_cat.head()
Out[41]:
In [42]:
#check number of NAs
health_cat.isnull().sum().sum()
Out[42]:
In [43]:
health_cat.shape
Out[43]:
In [44]:
#create reduced numerical dataframe
health_num = health_na[['YRSWRKPA', 'ASISLEEP', 'AHEIGHT', 'AWEIGHTP', 'BMI', 'BEDDAYR', 'CLCKTP', 'AGE_P', 'AHCNOYR2', 'LOCALL1B']]
In [45]:
health_num.head()
Out[45]:
In [46]:
health_num['LOCALL1B'].shape
Out[46]:
In [47]:
#check number of NAs
health_num.isnull().sum().sum()
Out[47]:
In [48]:
plt.show(health_num.boxplot(column=['YRSWRKPA', 'ASISLEEP', 'AHEIGHT','AWEIGHTP']))
In [49]:
plt.show(health_num.boxplot(column=['BMI']))
In [50]:
plt.show(health_num.boxplot(column=['BEDDAYR', 'CLCKTP', 'AGE_P', 'AHCNOYR2', 'LOCALL1B']))
In [51]:
# # z-score standard
# zscore = lambda x: ((x - x.mean()) / x.std()) if (x.dtypes==np.float64 or x.dtypes==np.int64) else x
# health_num_std = health_num.copy()
# health_num_std.apply(zscore).head()
In [52]:
#min max normalization of numerical variables
x = health_num.values #returns a numpy array
min_max_scaler = preprocessing.MinMaxScaler()
x_scaled = min_max_scaler.fit_transform(x)
health_num_norm = pd.DataFrame(x_scaled, columns=health_num.columns)
health_num_norm.head()
Out[52]:
In [53]:
health_num_norm['LOCALL1B'].shape
Out[53]:
In [54]:
#check number of NAs
health_num_norm.isnull().sum()
Out[54]:
In [55]:
plt.show(health_num_norm.boxplot(column=['YRSWRKPA', 'ASISLEEP', 'AHEIGHT', 'AWEIGHTP', 'BMI', 'BEDDAYR', 'CLCKTP', 'AGE_P', 'AHCNOYR2', 'LOCALL1B']))
In [56]:
#check for correlatin between numerical variables
health_num_norm.corr()
Out[56]:
In [57]:
pd.plotting.scatter_matrix(health_num_norm, figsize=(16, 16), alpha = 0.2)
plt.show()
In [58]:
plt.figure(figsize=(10,10))
sns.heatmap(health_num_norm.corr(), cmap='BuGn')
Out[58]:
In [59]:
#avr weight and avr height are higly correlated
#age and age on the job is correlated
#drop avr weight and age
health_num_norm_drop = health_num_norm.drop("AWEIGHTP",axis=1)
In [60]:
health_num_norm_drop['LOCALL1B'].shape
Out[60]:
In [61]:
#check number of NAs
health_num_norm_drop.isnull().sum().sum()
Out[61]:
In [62]:
#drop age
health_num_norm_drop = health_num_norm_drop.drop("AGE_P",axis=1)
In [63]:
#check number of NAs
health_num_norm_drop.isnull().sum().sum()
Out[63]:
In [64]:
#Plot correlation
sns.heatmap(health_num_norm_drop.corr(), cmap='BuGn')
Out[64]:
In [65]:
#re-run correlations
health_num_norm_drop.corr()
Out[65]:
In [66]:
#save health_num_norm_drop to csv file
health_num_norm_drop.to_csv('health_num_norm_drop.csv', sep=",")
In [67]:
health_num_norm_drop.shape
Out[67]:
In [68]:
#force conversion to categorical varibles
health_cat = health_cat.astype('category')
In [69]:
print(health_cat.dtypes.to_string())
In [70]:
health_cat['ALCSTAT'].shape
Out[70]:
In [71]:
health_dummy = pd.get_dummies(health_cat,drop_first=True)
health_dummy.reset_index(drop=True, inplace=True)
health_dummy.head()
Out[71]:
In [72]:
health_dummy['ALCSTAT_10'].shape
Out[72]:
In [73]:
#check number of NAs
health_dummy.isnull().sum().sum()
Out[73]:
In [74]:
columnNamesDummy = list(health_dummy.head(0))
print(columnNamesDummy)
In [75]:
#save dummy variables to csv file
health_dummy.to_csv('health_dummy.csv', sep=",")
In [76]:
health_dummy.shape
Out[76]:
In [77]:
health_num.shape
Out[77]:
In [78]:
health_num.isnull().sum().sum()
Out[78]:
In [79]:
#reset index for health_num
health_num.reset_index(drop=True, inplace=True)
In [80]:
#merge categorical dummy with numerical non-normalized
health_cleaned=pd.concat([health_dummy,health_num],axis=1,join_axes=[health_dummy.index])
health_cleaned.tail()
Out[80]:
In [81]:
# #merge categorical dummy with numerical normalized
# health_cleaned=pd.concat([health_dummy,health_num_norm_drop],axis=1,join_axes=[health_dummy.index])
# health_cleaned.tail()
In [82]:
health_cleaned['LOCALL1B'].shape
Out[82]:
In [83]:
health_cleaned['SEX_2'].shape
Out[83]:
In [84]:
health_cleaned.shape
Out[84]:
In [85]:
#save health cleaned to csv file
health_cleaned.to_csv('health_cleaned.csv', sep=",")
In [86]:
#Y ->'DBHVPAY_2.0' = Told to increase physical activity, past 12 m
In [87]:
#create dataframe
health_cl = health_cleaned[['DBHVPAY_2.0', 'SEX_2', 'R_MARITL_2', 'R_MARITL_3', 'R_MARITL_4', 'MRACRPI2_2', 'MRACRPI2_3', 'MRACRPI2_4', 'REGION_2', 'REGION_3', 'REGION_4', 'PAR_STAT_2', 'PAR_STAT_3', 'DOINGLWA_2.0', 'DOINGLWA_3.0', 'DOINGLWA_4.0', 'DOINGLWA_5.0', 'SUPERVIS_2.0', 'WRKCATA_2.0', 'WRKCATA_3.0', 'WRKCATA_4.0', 'WRKCATA_5.0', 'WRKCATA_6.0', 'HOURPDA_2.0', 'PDSICKA_2.0', 'WRKLYR4_1.0', 'WRKLYR4_2.0', 'HYPEV_2.0', 'HYBPLEV_2.0', 'HYBPLEV_3.0', 'HYBPLEV_4.0', 'HYBPLEV_5.0', 'CHLEV_2.0', 'CHDEV_2.0', 'MIEV_2.0', 'STREV_2.0', 'COPDEV_2.0', 'AASMEV_2.0', 'ULCEV_2.0', 'CANEV_2.0', 'DBHVCLY_2.0', 'DBHVWLY_2.0', 'DBHVPAN_2.0', 'DBHVCLN_2.0', 'DBHVWLN_2.0', 'DIBREL_2.0', 'DIBEV1_3.0', 'DIBPRE2_2.0', 'EPILEP1_2.0', 'AHAYFYR_2.0', 'SINYR_2.0', 'CBRCHYR_2.0', 'KIDWKYR_2.0', 'LIVYR_2.0', 'JNTSYMP_2.0', 'ARTH1_2.0', 'PAINECK_2.0', 'PAINLB_2.0', 'PAINFACE_2.0', 'AMIGR_2.0', 'ACOLD2W_2.0', 'AINTIL2W_2.0', 'AHEARST1_2.0', 'AHEARST1_3.0', 'AHEARST1_4.0', 'AHEARST1_5.0', 'AHEARST1_6.0', 'AVISION_2.0', 'VIM_GLEV_2.0', 'VIM_MDEV_2.0', 'VIMGLASS_2.0', 'AVISACT_2.0', 'CHPAIN6M_2.0', 'CHPAIN6M_3.0', 'CHPAIN6M_4.0', 'AHSTATYR_2.0', 'AHSTATYR_3.0', 'FLA1AR_2', 'FLA1AR_3', 'SPECEQ_2.0', 'ALC1YR_2.0', 'CIGAREV2_2.0', 'ECIGEV2_2.0', 'SMKSTAT2_2.0', 'SMKSTAT2_3.0', 'SMKSTAT2_4.0', 'APLKIND_2.0', 'APLKIND_3.0', 'APLKIND_4.0', 'APLKIND_5.0', 'APLKIND_6.0', 'AWORPAY_2.0', 'AWORPAY_3.0', 'ADNLONG2_1.0', 'ADNLONG2_2.0', 'ADNLONG2_3.0', 'ADNLONG2_4.0', 'ADNLONG2_5.0', 'ASRGYR_2.0', 'AMDLONGR_1.0', 'AMDLONGR_2.0', 'AMDLONGR_3.0', 'AMDLONGR_4.0', 'AMDLONGR_5.0', 'HIT1A_2.0', 'HIT2A_2.0', 'HIT3A_2.0', 'HIT4A_2.0', 'FLUVACYR_2.0', 'LIVEV_2.0', 'ASICPUSE_2.0', 'ASICPUSE_3.0', 'ASICPUSE_4.0', 'ASIRETR_2.0', 'ASIRETR_3.0', 'ASIRETR_4.0', 'ASIMEDC_2.0', 'ASIMEDC_3.0', 'ASIMEDC_4.0', 'ASISTLV_2.0', 'ASISTLV_3.0', 'ASISTLV_4.0', 'ASICNHC_2.0', 'ASICNHC_3.0', 'ASICNHC_4.0', 'AWEBUSE_2.0', 'YTQU_YG1_2.0', 'ALCSTAT_2', 'ALCSTAT_3', 'ALCSTAT_5', 'ALCSTAT_6', 'ALCSTAT_7', 'ALCSTAT_8', 'ALCSTAT_9', 'ALCSTAT_10', 'YRSWRKPA', 'ASISLEEP', 'AHEIGHT', 'BMI', 'BEDDAYR', 'CLCKTP', 'AHCNOYR2', 'LOCALL1B']]
health_cl['DBHVPAY_2.0'].head()
Out[87]:
In [88]:
#recode 'DBHVPAY_2.0' so 1 means Yes
physical = pd.get_dummies(health_cl['DBHVPAY_2.0'])
physical.columns = ["phys_yes", "phys_no"]
physical = physical.drop(columns=['phys_no'])
physical.head()
Out[88]:
In [89]:
#repalce 'phys_yes' with 'DBHVPAY_YES'
health_tree=pd.concat([physical,health_cl],axis=1,join_axes=[physical.index])
health_tree = health_tree.drop(columns=['DBHVPAY_2.0'])
health_tree.columns.values[0] = "DBHVPAY_YES"
health_tree.head()
Out[89]:
In [90]:
health_tree_drop = health_tree
In [91]:
#column names
health_names = health_tree_drop.columns.values
health_names
Out[91]:
In [92]:
y = health_tree_drop['DBHVPAY_YES'] # y variable
X = health_tree_drop[health_names[1:]]
X.head()
Out[92]:
In [93]:
y.head
Out[93]:
In [94]:
#create train and test sets
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=99874)
print "X_train", X_train.shape
print "y_train", y_train.shape
print "X_test", X_test.shape
print "y_test", y_test.shape
In [95]:
############################################################################################
#Decistion Tree
############################################################################################
In [96]:
treeclf = tree.DecisionTreeClassifier(criterion='entropy', min_samples_split=3)
treeclf = treeclf.fit(X_train, y_train)
In [97]:
#A versatile function to measure performance of a classification model
from sklearn import metrics
def measure_performance(X, y, clf, show_accuracy=True, show_classification_report=True, show_confussion_matrix=True):
y_pred = clf.predict(X)
if show_accuracy:
print "Accuracy:{0:.3f}".format(metrics.accuracy_score(y, y_pred)),"\n"
if show_classification_report:
print "Classification report"
print metrics.classification_report(y, y_pred),"\n"
if show_confussion_matrix:
print "Confussion matrix"
print metrics.confusion_matrix(y, y_pred),"\n"
In [98]:
measure_performance(X_test, y_test, treeclf, show_confussion_matrix=True, show_classification_report=True)
In [99]:
#predictions
treepreds_test = treeclf.predict(X_test)
print treepreds_test
In [100]:
#classification report
print(classification_report(y_test, treepreds_test))
In [101]:
treecm = confusion_matrix(y_test, treepreds_test)
print treecm
In [102]:
print treeclf.score(X_test, y_test)
print treeclf.score(X_train, y_train)
In [103]:
export_graphviz(treeclf,out_file='tree.dot', feature_names=X_train.columns, class_names=["No","Yes"])
with open("tree.dot") as f:
dot_graph = f.read()
graphviz.Source(dot_graph)
Out[103]:
In [104]:
system(dot -Tpng tree.dot -o dtree.png)
Out[104]:
In [105]:
from IPython.display import Image
Image(filename='dtree.png', width=800)
Out[105]:
In [106]:
########################################################################################
#Feature Selection
########################################################################################
In [107]:
#Select the top 30% of the most important features, using a chi2 test
fs = feature_selection.SelectPercentile(feature_selection.chi2, percentile=30)
X_train_fs = fs.fit_transform(X_train, y_train)
In [108]:
np.set_printoptions(suppress=True, precision=2, linewidth=320)
print fs.get_support()
print fs.scores_
In [109]:
X.columns.values
Out[109]:
In [110]:
col_fs = X.columns[fs.get_support()].values
print X.columns[fs.get_support()].values
In [111]:
for i in range(len(X.columns.values)):
if fs.get_support()[i]:
print X.columns.values[i],'\t', fs.scores_[i]
In [112]:
print X_train_fs
In [113]:
dt = tree.DecisionTreeClassifier(criterion='entropy')
dt.fit(X_train_fs, y_train)
X_test_fs = fs.transform(X_test)
measure_performance(X_test_fs, y_test, dt, show_confussion_matrix=False, show_classification_report=True)
In [114]:
#export_graphviz(dt,out_file='tree2.dot', feature_names=X_train_fs.columns, class_names=["No","Yes"])
export_graphviz(dt,out_file='tree2.dot', feature_names=col_fs, class_names=["No","Yes"])
with open("tree2.dot") as f:
dot_graph = f.read()
graphviz.Source(dot_graph)
Out[114]:
In [115]:
system(dot -Tpng tree2.dot -o dtree2.png)
Out[115]:
In [116]:
from IPython.display import Image
Image(filename='dtree2.png', width=800)
Out[116]:
In [117]:
#Find the best percentile for feature selection using cross-validation
from sklearn import cross_validation
dt = tree.DecisionTreeClassifier(criterion='entropy')
percentiles = range(1, 100, 5)
results = []
for i in range(1, 100, 5):
fs = feature_selection.SelectPercentile(feature_selection.chi2, percentile=i)
X_train_fs = fs.fit_transform(X_train, y_train)
scores = cross_validation.cross_val_score(dt, X_train_fs, y_train, cv=5)
print i,scores.mean()
results = np.append(results, scores.mean())
# optimal_percentile = np.where(results == results.max())[0]
# print "Optimal percentile of features:{0}".format(percentiles[optimal_percentile]), "\n"
# optimal_num_features = int(percentiles[optimal_percentile]*len(X.columns)/100)
# print "Optimal number of features:{0}".format(optimal_num_features), "\n"
optimal_percentil = int(np.where(results == results.max())[0])
print "Optimal percentile of features is: ", percentiles[optimal_percentil]
#optimal_num_features = percentiles[optimal_percentile]*len(X.columns)/100
#optimal_num_features = int(math.floor(percentiles[optimal_percentil]*X.shape[1]/100))
optimal_num_features = int((percentiles[optimal_percentil]*X.shape[1]/100))
print "Optimal number of features is: ", optimal_num_features
# Plot percentile of features VS. cross-validation scores
import pylab as pl
pl.figure()
pl.xlabel("Percentage of features selected")
pl.ylabel("Cross validation accuracy")
pl.plot(percentiles,results)
Out[117]:
In [118]:
#Impact of max-depth
In [119]:
from sklearn.cross_validation import KFold
def calc_params(X, y, clf, param_values, param_name, K):
# Convert input to Numpy arrays
X = np.array(X)
y = np.array(y)
# initialize training and testing scores with zeros
train_scores = np.zeros(len(param_values))
test_scores = np.zeros(len(param_values))
# iterate over the different parameter values
for i, param_value in enumerate(param_values):
print param_name, ' = ', param_value
# set classifier parameters
clf.set_params(**{param_name:param_value})
# initialize the K scores obtained for each fold
k_train_scores = np.zeros(K)
k_test_scores = np.zeros(K)
# create KFold cross validation
cv = KFold(len(X), K, shuffle=True, random_state=0)
# iterate over the K folds
for j, (train, test) in enumerate(cv):
# fit the classifier in the corresponding fold
# and obtain the corresponding accuracy scores on train and test sets
clf.fit([X[k] for k in train], y[train])
k_train_scores[j] = clf.score([X[k] for k in train], y[train])
k_test_scores[j] = clf.score([X[k] for k in test], y[test])
# store the mean of the K fold scores
train_scores[i] = np.mean(k_train_scores)
test_scores[i] = np.mean(k_test_scores)
# plot the training and testing scores in a log scale
plt.plot(param_values, train_scores, label='Train', alpha=0.4, lw=2, c='b')
plt.plot(param_values, test_scores, label='X-Val', alpha=0.4, lw=2, c='g')
plt.legend(loc=7)
plt.xlabel(param_name + " values")
plt.ylabel("Mean cross validation accuracy")
# return the training and testing scores on each parameter value
return train_scores, test_scores
In [120]:
# Let's create an evenly spaced range of numbers in a specified interval
md = np.linspace(1, 40, 20)
md = np.array([int(e) for e in md])
print md
In [121]:
train_scores, test_scores = calc_params(X_train, y_train, dt, md, 'max_depth', 5)
In [122]:
#max_depth = 4 looks the best
In [123]:
#min number of samples allowed at a leaf node
In [124]:
msl = np.linspace(1, 30, 15)
msl = np.array([int(e) for e in msl])
dt = tree.DecisionTreeClassifier(criterion='entropy')
train_scores, test_scores = calc_params(X_train, y_train, dt, msl, 'min_samples_leaf', 5)
In [125]:
#min_samples_leaf should be between 10 and 15
In [126]:
#Grid Search allows us to more systemiatically explore different combinations of multiple parameters
In [127]:
#from sklearn.grid_search import GridSearchCV
from sklearn.model_selection import GridSearchCV
dt = tree.DecisionTreeClassifier()
# parameters = {
# 'criterion': ['entropy','gini'],
# 'max_depth': np.linspace(1, 20, 10),
# 'min_samples_leaf': np.linspace(10, 15, 5),
# 'min_samples_split': np.linspace(10, 15, 5)
# }
parameters = {
'criterion': ['entropy','gini'],
'max_depth': [1,2,3,4,5,6,8,10],
'min_samples_leaf': [2,4,6,8,10,12,14,16,18,20,],
'min_samples_split': [2,4,6,8,10,12,14,16,18,20]
}
gs = GridSearchCV(estimator=dt, param_grid=parameters, verbose=1, cv=5,refit=True, n_jobs = -1)
In [128]:
%time
gs.fit(X_train, y_train)
gs.best_params_, gs.best_score_
Out[128]:
In [129]:
# Note: best parameters results change most of the time
# [Parallel(n_jobs=1)]: Done 8000 out of 8000 | elapsed: 2.9min finished
# ({'criterion': 'gini',
# 'max_depth': 5,
# 'min_samples_leaf': 20,
# 'min_samples_split': 2},
# 0.8076197957580519)
dt = tree.DecisionTreeClassifier(criterion='gini', max_depth=5, min_samples_leaf=20, min_samples_split=2)
dt.fit(X_train, y_train)
measure_performance(X_test, y_test, dt, show_confussion_matrix=False, show_classification_report=True)
In [130]:
from sklearn.tree import export_graphviz
export_graphviz(dt,out_file='treefinal.dot', feature_names=X_train.columns)
In [131]:
import graphviz
with open("treefinal.dot") as f:
dot_graph = f.read()
graphviz.Source(dot_graph)
Out[131]:
In [132]:
system(dot -Tpng treefinal.dot -o treefinal.png)
Out[132]:
In [133]:
from IPython.display import Image
Image(filename='treefinal.png', width=800)
Out[133]:
In [134]:
print(X_train.columns.values)
In [135]:
dt = tree.DecisionTreeClassifier(criterion='entropy', max_depth=4, min_samples_leaf=10, min_samples_split=2)
dt.fit(X_train, y_train)
measure_performance(X_test, y_test, dt, show_confussion_matrix=True, show_classification_report=True)
In [136]:
dt
Out[136]:
In [137]:
dt.get_params
Out[137]:
In [138]:
from sklearn.tree import export_graphviz
export_graphviz(dt,out_file='treefinal2.dot', feature_names=X_train.columns)
In [139]:
import graphviz
with open("treefinal2.dot") as f:
dot_graph = f.read()
graphviz.Source(dot_graph)
Out[139]:
In [140]:
system(dot -Tpng treefinal2.dot -o treefinal2.png)
Out[140]:
In [141]:
from IPython.display import Image
Image(filename='treefinal2.png', width=800)
Out[141]:
In [142]:
#variables returned by decision tree
#'DBHVCLY_2.0','BMI','DBHVPAN_2.0','DBHVWLY_2.0','CANEV_2.0','ASISTLV_4.0','AHCNOYR2','BEDDAYR',
#'DBHVCLN_2.0','AMDLONGR_1.0','AHCNOYR2','FLA1AR_2'
In [143]:
#create new data frame
#health_cleaned
health_tree.tail()
Out[143]:
In [144]:
health_tree= pd.DataFrame(health_tree)
In [145]:
health_tree.head(1)
Out[145]:
In [146]:
hn = health_tree[['DBHVPAY_YES','DBHVCLY_2.0','BMI','DBHVPAN_2.0','DBHVWLY_2.0','CANEV_2.0','ASISTLV_4.0','AHCNOYR2','BEDDAYR','DBHVCLN_2.0','AMDLONGR_1.0','AHCNOYR2','FLA1AR_2']]
hn.head()
Out[146]:
In [147]:
y = hn['DBHVPAY_YES'] # y variable
X = hn[['DBHVCLY_2.0','BMI','DBHVPAN_2.0','DBHVWLY_2.0','CANEV_2.0','ASISTLV_4.0','AHCNOYR2','BEDDAYR','DBHVCLN_2.0','AMDLONGR_1.0','AHCNOYR2','FLA1AR_2']]
X.head()
Out[147]:
In [148]:
#create train and test sets
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=9874)
print "X_train", X_train.shape
print "y_train", y_train.shape
print "X_test", X_test.shape
print "y_test", y_test.shape
In [149]:
##############################################################################################################
#Logistic Regression
##############################################################################################################
In [150]:
from sklearn.linear_model import LogisticRegression
In [151]:
#logistic regression with recursive feature elimination using cross validation
LogReg = LogisticRegression()
refcv = RFECV(LogReg, cv = 10, step = 1, n_jobs = -1)
fitcv = refcv.fit(X_train, y_train)
In [152]:
y_pred = fitcv.predict(X_test)
y_pred
Out[152]:
In [153]:
from sklearn.metrics import confusion_matrix
confusion_matrix = confusion_matrix(y_test, y_pred)
confusion_matrix
Out[153]:
In [154]:
print(classification_report(y_test, y_pred))
In [155]:
#np.set_printoptions(threshold=np.inf)
np.set_printoptions(threshold=1000)
print "", y_pred
In [156]:
print("Accuracy:",metrics.accuracy_score(y_test, y_pred))
print("Precision:",metrics.precision_score(y_test, y_pred))
print("Recall:",metrics.recall_score(y_test, y_pred))
In [157]:
y_pred_proba = fitcv.predict_proba(X_test)[::,1]
fpr, tpr, _ = metrics.roc_curve(y_test, y_pred_proba)
auc = metrics.roc_auc_score(y_test, y_pred_proba)
plt.plot(fpr,tpr,label="auc="+str(auc))
plt.legend(loc=4)
plt.show()
print "Area Under Curve AUC is ", auc
In [158]:
fitcv.score(X_train, y_train)
Out[158]:
In [159]:
fitcv.score(X_test, y_test)
Out[159]:
In [160]:
print("Num Features: %d") % fitcv.n_features_
In [161]:
print("Selected Features: %s") % fitcv.support_
In [162]:
lr_selected = fitcv.support_
lr_selected
Out[162]:
In [163]:
lr_rank = fitcv.ranking_
print("Feature Ranking: %s") % fitcv.ranking_
In [165]:
lr_cols = X_train.columns.values
lr_cols
Out[165]:
In [167]:
for i in range(len(lr_cols)):
#if lr_selected[i]:
print lr_cols[i],'\t', lr_rank[i]
In [168]:
#############################################################################################################
#logistic regression with cross validation using recursive feature elimination
#############################################################################################################
In [169]:
from sklearn.linear_model import LogisticRegressionCV
from sklearn.feature_selection import RFE
In [170]:
X_train.tail()
Out[170]:
In [171]:
LogRegCV = LogisticRegressionCV(cv=10, random_state=123231, multi_class='ovr')
In [172]:
rfe = RFE(LogRegCV,5)
In [173]:
fit = rfe.fit(X_train, y_train)
In [174]:
y_pred_cv = fit.predict(X_test)
y_pred_cv
Out[174]:
In [175]:
from sklearn.metrics import confusion_matrix
confusion_matrix = confusion_matrix(y_test, y_pred_cv)
confusion_matrix
Out[175]:
In [176]:
lg = classification_report(y_test, y_pred_cv)
print lg
#print(classification_report(y_test, y_pred_cv))
In [177]:
import pylab as plt
%matplotlib inline
plt.matshow(confusion_matrix)
plt.title('Confusion matrix')
plt.colorbar()
plt.ylabel('Actual')
plt.xlabel('Predicted')
plt.show()
In [178]:
print("Accuracy:",metrics.accuracy_score(y_test, y_pred_cv))
print("Precision:",metrics.precision_score(y_test, y_pred_cv))
print("Recall:",metrics.recall_score(y_test, y_pred_cv))
In [179]:
y_pred_proba_cv = fit.predict_proba(X_test)[::,1]
fpr, tpr, _ = metrics.roc_curve(y_test, y_pred_proba_cv)
auc = metrics.roc_auc_score(y_test, y_pred_proba_cv)
plt.plot(fpr,tpr,label="auc="+str(auc))
plt.legend(loc=4)
plt.show()
print "Area Under Curve AUC is ", auc
In [180]:
fit.score(X_train, y_train)
Out[180]:
In [181]:
fit.score(X_test, y_test)
Out[181]:
In [182]:
print("Num Features: %d") % fit.n_features_
In [183]:
print("Selected Features: %s") % fit.support_
In [184]:
lr_selected = fit.support_
lr_selected
Out[184]:
In [185]:
lr_rank = fit.ranking_
print("Feature Ranking: %s") % fit.ranking_
In [186]:
lr_cols = X_train.columns.values
lr_cols
Out[186]:
In [187]:
for i in range(len(lr_cols)):
#if lr_selected[i]:
print lr_cols[i],'\t', lr_rank[i]
In [188]:
# varaibles selected: DBHVCLY_2.0, DBHVPAN_2.0, DBHVWLY_2.0, AMDLONGR_1.0, FLA1AR_2
In [189]:
###############################################################################################
# Logistic Regression using stasmodels.org api
#################################################################################################
In [190]:
import statsmodels.api as sm
In [191]:
hn.head()
Out[191]:
In [192]:
ys = health_tree_drop['DBHVPAY_YES'] # y variable
Xs = health_tree_drop[health_names[1:]]
In [193]:
ys.tail()
Out[193]:
In [194]:
Xs.dtypes
Out[194]:
In [195]:
Xs.dtypes
Out[195]:
In [196]:
logit_model=sm.Logit(ys,Xs)
result=logit_model.fit()
print(result.summary2())
In [197]:
cols=['SEX_2', 'MRACRPI2_4', 'WRKCATA_2.0', 'WRKCATA_5.0', 'HOURPDA_2.0', 'PDSICKA_2.0', 'HYPEV_2.0', 'HYBPLEV_5.0', 'CHLEV_2.0', 'CANEV_2.0', 'DBHVCLY_2.0', 'DBHVWLY_2.0', 'DBHVPAN_2.0', 'DBHVCLN_2.0', 'DBHVWLN_2.0', 'DIBREL_2.0', 'DIBPRE2_2.0', 'AHEARST1_2.0', 'AHEARST1_4.0', 'AHEARST1_6.0', 'VIMGLASS_2.0', 'AVISACT_2.0', 'FLA1AR_2', 'APLKIND_5.0', 'AMDLONGR_4.0', 'AMDLONGR_5.0', 'HIT1A_2.0', 'HIT4A_2.0', 'ASICPUSE_2.0', 'ASICPUSE_4.0', 'ASIMEDC_2.0', 'ASISTLV_4.0', 'AWEBUSE_2.0', 'YTQU_YG1_2.0', 'ASISLEEP', 'BMI', 'AHCNOYR2']
cols[0:5]
Out[197]:
In [198]:
Xs2= health_tree_drop[cols]
Xs2[0:5]
Out[198]:
In [199]:
ys2= ys
ys2.head()
Out[199]:
In [200]:
logit_model2=sm.Logit(ys2,Xs2)
result2=logit_model2.fit()
print(result2.summary2())
In [201]:
#crate test train
X_train2, X_test2, y_train2, y_test2 = train_test_split(Xs2, ys2, test_size=0.2, random_state=2323)
LogReg3 = LogisticRegression()
LogReg3.fit(X_train2, y_train2)
Out[201]:
In [202]:
#predict y
y_pred2 = LogReg3.predict(X_test2)
print('Accuracy : {:.2f}'.format(LogReg3.score(X_test2, y_test2)))
In [203]:
#conf matrix
from sklearn.metrics import confusion_matrix
confusion_matrix3 = confusion_matrix(y_test2, y_pred2)
print(confusion_matrix3)
In [204]:
lg2 = classification_report(y_test2, y_pred2)
print lg2
In [205]:
y_pred_proba2 = LogReg3.predict_proba(X_test2)[::,1]
fpr, tpr, _ = metrics.roc_curve(y_test2, y_pred_proba2)
auc = metrics.roc_auc_score(y_test2, y_pred_proba2)
plt.plot(fpr,tpr,label="auc="+str(auc))
plt.legend(loc=4)
plt.show()
print "Area Under Curve AUC is ", auc
In [206]:
##############################################################################################################
#LDA
##############################################################################################################
In [207]:
ldclf = LinearDiscriminantAnalysis()
ldclf = ldclf.fit(X_train, y_train)
print "Score on Training: ", ldclf.score(X_train, y_train)
print "Score on Test: ", ldclf.score(X_test, y_test)
In [208]:
# 10-fold cross validation
cv_scores = cross_validation.cross_val_score(ldclf, X, y, cv=5)
cv_scores
Out[208]:
In [209]:
print("Overall Accuracy: %0.2f (+/- %0.2f)" % (cv_scores.mean(), cv_scores.std() * 2))
In [222]:
##############################################################################################################
#Clustering numerical variables
##############################################################################################################
In [223]:
health_cl = health_tree
print health_cl.shape
health_cl.tail()
Out[223]:
In [224]:
health_cl_num = health_cl[['YRSWRKPA','ASISLEEP','AHEIGHT','BMI','BEDDAYR','CLCKTP','AHCNOYR2','LOCALL1B']]
health_cl_num.tail()
Out[224]:
In [225]:
#min max normalization of numerical variables
x = health_cl_num.values #returns a numpy array
min_max_scaler = preprocessing.MinMaxScaler()
x_scaled = min_max_scaler.fit_transform(x)
health_cl_num_norm = pd.DataFrame(x_scaled, columns=health_cl_num.columns)
health_cl_num_norm.columns = ['YRSWRKPA_NORM','ASISLEEP_NORM','AHEIGHT_NORM','BMI_NORM','BEDDAYR_NORM','CLCKTP_NORM','AHCNOYR2_NORM','LOCALL1B_NORM']
health_cl_num_norm.head()
Out[225]:
In [226]:
health_cl_num_norm.dtypes
Out[226]:
In [227]:
#Xcl = np.array(health_cl2.drop(['DBHVPAY_YES'], 1).astype(float64))
#Xcl = np.array(health_cl2.drop(['DBHVPAY_YES'], 1))
Xcl = np.array(health_cl_num_norm)
In [228]:
Xcl[0:5]
Out[228]:
In [231]:
#
kmeans = KMeans(n_clusters=5, max_iter=500, verbose = 0)
kmeans.fit(Xcl)
Out[231]:
In [233]:
#categorical dataset
health_cl1 = health_cl.iloc[:,0:135]
health_cl1.tail()
Out[233]:
In [234]:
#merge
health_cl2 = pd.concat([health_cl1,health_cl_num_norm],axis=1,join_axes=[health_cl1.index])
health_cl2.tail()
Out[234]:
In [235]:
#Xcl = np.array(health_cl2.drop(['DBHVPAY_YES'], 1).astype(float64))
#Xcl = np.array(health_cl2.drop(['DBHVPAY_YES'], 1))
Xcl = np.array(health_cl_num_norm)
In [236]:
Xcl[0:5]
Out[236]:
In [237]:
ycl = np.array(health_cl2['DBHVPAY_YES'])
In [238]:
ycl[0:5]
Out[238]:
In [239]:
#
kmeans = KMeans(n_clusters=5, max_iter=500, verbose = 1)
kmeans.fit(Xcl)
Out[239]:
In [240]:
clusters = kmeans.predict(Xcl)
print clusters
In [241]:
clusters.shape
Out[241]:
In [242]:
ycl.shape
Out[242]:
In [243]:
print completeness_score(ycl,clusters)
In [244]:
print homogeneity_score(ycl,clusters)
In [ ]:
#this clustering model cannot be used, it was an experiment for clustering numerical data with binary y variable.